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1 Introduction 

In this article we present a novel way to estimate the amounts of traffic on the Origin- 
Destination couples (OD couples). This new approach combines together a routing 
algorithm based on the principle of the shortest path and a recent technique of stochas- 
tic optimization called Cross-Entropy. The CE method was built at the origin, to 
tackle problems of rare-event simulation. However, its inventor, R. Rubinstein, re- 
alized soon that the underlying idea should be applied efficiently to combinatorial and 
multi-extremal optimization problems. 

In a final part, we adapt a particular filtering algorithm in order to be able to dynami- 
cally estimate the evolution of the traffic on the OD couples. 

The aim of this report is to highlight rather original ideas, however the choices of the 
prior distributions and some specific parameters may be quite arbitrary. 

2 A brief presentation of the CE method 

2.1 Rare-Event Simulation 

Let X — (Xi, Xn), be a random vector taking values in a space called X. Let 
{/(.; v)} be a family of parametric densities defined on the space X, with respect to 
the Lebesgue measure. For any measurable function H, we can define: 



E[H{X)] = j H(x) f(x;v) dx 



The performance function will be called S : X — > R. For a fixed level 7, we are 
interessed in the probabilty of the event defined below: 

I = Pu(S(X) > 7) = E M [1 {S(X) > 7} ] . 
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If this probability is very small, for example not more than 10 -5 , the set {S(x) > 
7} will be called a rare event. 

A straightforward way to estimate I may be to use crude Monte-Carlo simulation: let 
(Xi, .... Xn) be a sample drawn from the density /(.; u). Then, the estimator, 



i=l 



N ^ 1 {S(x)> 1 } 



is an unbiased estimator of I. However, if {S(x) > 7} is a rare event, many 
indicator functions will remain equal to zero. As a result, we will be forced to simulate 
huge samples, which is rather costly and difficult to put in application. 
Another way to get an estimate of I might be to use importance sampling. We should 
draw (Xi, Xn) from a density g , defined on the space X. This density is nothing 
else than a mere change of measure. The estimator then becomes: 

1 - 1 V 1 ffiig) m 

1 - 1 ^™>7}-^f- (D 

The optimal density g, is defined by: 

g*(x) = 1{S( ^] } fiX;U) . (2) 
Substituting (0 in Q, we get: 

But, I is a constant. As a result, the estimator defined in Q has zero variance. Nev- 
ertheless, g* depends on the unknown parameter /. The idea is in fact, to choose g in a 
parametric family of densities {/(.;«)}. The problem is now to determine the optimal 
parameter v, such that the distance between g* and /(.; v) should be minimized. 
A well-known "distance" between two densities g and h, is the Kullback-Leibler "dis- 
tance": 



T>(g,h) = E g [ln ||^|] = / .</(•'■)/»//(•'') - / uU) In In.r) ,l.r . (3) 



Minimizing the Kullback-Leibler distance between g* and /(.; v), is equivalent to 
solve the following problem: 



argmax / g*(x) In f(x;v) dx . (4) 
" J 

Substituting (01 in Q, we get: 

argmaxD(u) = argmax ~E U [1 {S ( X ) < 7 } In f(X; v)] . (5) 
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But, in fact we can estimate v*, by solving the following stochastic program: 

1 N 

argmax D(v) = argmax — 2J [ 1 {S(x i ) > 7 } Jn f(Xi;v)] . (6) 

i=l 

If D is convex and differentiable in v, we just have to solve the following problem: 

1 N 

]V E ^Wxo > 7} V, /(X i5 «)] = . (7) 

i=l 

The solution can often be calculated analytically, which is one of the great advan- 
tages of this approach. 

2.2 Application of the CE to optimization 

Usually, in the field of optimization we try to solve problems of the form: 

S(x*) = 7* = argmax S(x) . (8) 

The genius of the Cross-Entropy method lies in the fact that it is possible to asso- 
ciate with each optimization problem of the form (JgJ, a problem of estimation, called 
associated stochastic problem (ASP). We will start by defining a collection of indicator 
functions {l{s(x) < -y}}y e K, on the space X. Then, we will define a parametric fam- 
ily of densities {/(.; v), v € V} on the space X. Let u £ V. We will associate with 
JHJ, the following stochastic estimation problem: 



1(7) - P U (S(X) > 7 ) = E 1 {s(x)>7> /(*;«) = E M [l { s(x)> 7} ] , (9) 

If 7 = 7*, a natural estimator of the reference parameter v* is: 

1 ^ 

u* = argmax — ^ > 7 } In f(Xf, v) , (10) 

where the Xj are drawn from the density f(.,u). If 7 is very close to 7*, then 
/(.;«*) assigns most of its probability mass close to x* . In fact, in this case, we will 
have to choose u so that P U (S(X) > 7) is not too small. We can infer that u and 7 
are closely linked. 

We will use a two level procedure. Indeed, we will construct two sequences 71, ...,7t 
and ubj^i, ■■■,vt such that jx is close to 7* and v"t is such that the density assigns 
most of its mass in the state which maximizes the performance. 
The algorithm follows a two-step strategy: 

Algorithm 1. • Define v = u. Set t = 1. 
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• Generate a sample (X\, Xn) ~ /(.; v t -\). Compute the (1 — p)-quantile 
of the performance S, which can be estimated by: 

it = <%i-p);v] ■ 

• Use the same sample (Xi, Xn), to solve ©. Call the solution Vt- 

• If for some t > d, d fixed, 

it = ■■■ = it-d , 

Stop; otherwise set t = t + 1 and reiterate from Step 2. 

3 The Model 

The network we study is composed of p nodes and n arcs. At first, we will suppose that 
there exists an arc between each couple of nodes. Furthermore, we make the difference 
between the two couples and (j, i), V i,j 6 {1, ••-,«}, i ^ j.. Consequently, 

n = p 2 — p . 

Recall that an arc is a directed link. In the rest of this article, we will make the hypoth- 
esis that the network is directed. 

Our work could be separated into two different parts. Firstly, we have to deal with 
a simulation part. In this section, we will choose an initial vector of the amounts of 
traffic on the OD couples. Our aim will be to minimize the global sum of the costs, 
which are associated to each arc of the network. Secondly, we will have to cope with 
an estimation part. Indeed, we will suppose that the initial costs are equal to those 
obtained by the simulation part. The idea then, will be to find the estimator X(t) from 
which we could infer an arc estimator Y(t) minimizing the distance to the vector Y(t), 
obtained in the simulation part. This raises the crucial question of the identifiability of 
the vector X(t). Is there unicity of the associated vector X(t) or, is it only an element 
in a vast variety? 

3.1 Simulation 

In this part, we associate a cost function to the network, which means that we give a cost 
to each arc of the network. This cost function is drawn from a parametric families of 
densities. Which means that we have to determine the optimal parameter of the density 
function. What's more, this cost may or, may not be, proportionnal to the amount of 
traffic on each arc. However, it is more unconventional to suppose that the costs depend 
on the arc traffic. Let Y(t) be the vector which contains the amount of traffic on each 
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arc of the network and C (t) the vector which represents the costs associated with each 
arc. In the most general case, we have: 

C(t) = F(Y(t)) , (11) 

where F is supposed to be continuous and differentiable. 



C(i) ~ f(-;vt) 



(+!)- 

t eration 




XO, UNDER SOME CONSTRAINTS 

■ L. -li. -i. iL. 

:■ . . . . . . ... . ....... . .... 

-> Y(t), C(t+1) = F(Y(t)) 



Figure 1: Simulation. 



Then, using these costs, we will use a routing algorithm based on the principle of 
the shortest path in order to find the shortest paths between each OD couple. A path is 
represented by the nodes or the arcs which it is made of. 

The principle of the routing algorithm we use is rather simple. The first observation to 
make is that every subpath of a shortest path is necessarily itself a shortest path. 
Let eij be the arc linking the nodes i andj. Ifthe path composed of the arcs {e^j, ej.fc, ek,i, e PtQ } 
is the shortest path linking node i to node q. Then, ejj must be the shortest path be- 
tween i and j. ej ; k must be the shortest path between j and k, and so on... 
ejj is called basic arc, iff it is the shortest path between i and j. Consequently, each 
shortest path must be composed exclusively of basic arcs. The aim of our routing algo- 
rithm is to substitute to each arc which is not basic, a basic arc. Let di j = C(t; i,j) 
be the distance or weight associated to the arc {i, j}, at the instant t. Let j be the indice 
of a node of the network, then: 
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\fi,k G Network — {j} , d^k <— min-fc^k, + dj^} ■ 

The algorithm tests all the couple of nodes (i, k), which are neighbors of j, while 
j takes each node of the network as its own value. 

Algorithm 2. Input: C(t), vector of the costs. 

• Ifdi t k > dij + dj^, do not change anything. 

• If di t k < di_j + dj k, create an arc linking i to k and associate the weight 

Output: the shortest paths between each OD couple. 

The algorithm also gives us the shortest distances associated to each OD couple. 
But, these distances are only rought estimators of the amounts of traffic on each OD 
couple. Indeed, more than one link, can be shared by different shortest paths linking 
different OD couples. As a result, the total amount of traffic generated by one OD 
couple usually represents only a fraction of the total traffic flowing through the arcs 
which composed the path. 



In order to solve this crucial problem, we will associate to the vector which contains 
the volumes of traffic on the OD couples, called Xo(t), an estimator of the volumes of 
traffic flowing through the arcs, which we will note Y(i). Indeed, if we use the routing 




Shortest Path 



OD couple ;J.5J 



Figure 2: Connection between OD couples and arcs. 
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algorithm, it is quite easy to deduce Y(t) from Xo(t). Our goal will be to solve the 
following optimization problem: 



in \Ci(t)\ . (12) 



mm 



3.2 Estimation 

In this part, the performance function is defined by: 



S(X(t)) = argmax ; . (13) 

x\t) \\Y(t)-Y(t)h 



1) -iteration 
(t+1 ) 



SIMULATED C(t) 



ROUTING 



ALGORITHM 



FTND AN ESTIMATOR OFX(t), SUCH THAT THE 
ASSOCIATED ESTIMATOR OF Y (t) MINIMIZES THE 
DISTANCE TO Y (t) 
fWFFR THE TRUNCATED COST 



Figure 3: Estimation 



The first observation is that S is an implicit function of X(t), e.g. we can't get any 
exact analytical expression of S. Consequently, we will have to resort to use simulation. 

What's more, we suppose that the vector X(t) is generated from an exponential 
density whose parameter is totally unknown. 

X{t) - £(A), A 6 M™ . (14) 

We apply the Cross-Entropy method to our problem. At each iteration, we generate 
a random sample (jW,..,!^) - £(A), A g R" . 

7 



Hypothesis: Each component of X(t), which represents the amount of traffic on 
an OD couple, will be supposed to be independent of the others. 



Practically, each random vector from the sample will be stocked in a big matrix. 



X(t) = 



(N) 



X. 



(1) 



X. 



(AT) 



X 



(1) 



X 



(JV) 



(15) 



The joined densities of the vectors X^\ j = 1, . .. , N, are typically of the form: 

n 



i=i 



= J] \i exp[ A, X™] 1 



?=1 



{min(Xl J> ) > 0} 



(16) 



As a result, we will have to solve the following problem: 



1 N 

v E l { s ( x«)> 7 - t} VWn/(XW;A) = 0. 



N 



(17) 



i=l 



After some computations, we get: 



J2iLi 1 {s(xm) > i t } 



A, = 



(i) 



, V j = 1, ...,n 



(18) 



Remark. 

If we generate the random vectors X i} i = 1, N from a truncated exponential, we 
can give some maximal boundaries on the OD volumes of traffic. 
Recall that a truncated exponential is of the form: 

,, . , , A exp[-A x] 
/(2;;A ' fe) = l-exp(-Afc) 1[ °- b]{x) - 

Under this assumption, we have to cope with the following system: 



N 



E 1 {S(XM) > it} X 

i=l - 

JV 

E 1 {S(X(*>) >it} 

i=l 



(i) 



+ 



Xj cxp(Aj bj)-l 



0, Vj e {l,...,n}. (19) 



This system is non-linear, that's why we use the well-known iterative Newton's 
method to solve it. 
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Figure 4: Simulation, 5 nodes network, 20 arcs. 



4 The problem of Identifiability 

We observe that the amount of traffic which flows through each arc of the network is 
equal to the sum of the amounts of traffic flowing on each OD couple which owns this 
arc in its shortest path. Remember that, thanks to the routing algorithm, we associate 
to each OD couple, a unique shortest path. Mathematically, we can express this remark 
under the following expression: 

n 

^ l{ej g OD couple number j'spath(t)} Xj(t) — Yi(t) , V i G {l,...,n}. (20) 

i=i 
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Figure 5: Simulation, 20 nodes network, 380 arcs. 



More generally, we get: 



Y(t) = A(X(t)) X(t) . (21) 

Indeed, in the most general case, the routing matrix A relies on the volumes of 
traffic flowing through each arc at the instant t. But, these arc volumes rely themselves 
on the OD volumes, X(t). As a result, the routing matrix A(t), is a function of X(t). 
A is uniquely made of binary elements: and 1. More explicitely, A(X(t);i, j) = 
1 iff, the arc numbered i belongs to the shortest path associated to the OD couple 
numbered j, at time t. 

What's more, the routing algorithm do not use every arc. Consequently, many rows of 
the routing matrix equal zero. The associated components in the arc vector Y (t) are at 
the same time, null. 

But, if we suppress the zero rows of A(X(t)) and the zero components of Y(t), this 
leads us to solve a rectangular system of equations. This system is under-determined, 
that's why we can't guarantee the existence of a unique solution. 



10 



We can conclude that there is no identifiability between the arc volumes Y(t), and 
the OD volumes X(t), at a given time. Indeed, if we take a fixed X(t), we get a 
unique associated Y(t), since the routing algorithm determines a unique shortest path 
between each OD couples. Reciprocally, if we take some fixed arc volumes, Y(t), we 
can't guarantee the unicity of the solutions of i|21|) . That's why, we can't assert that the 
associated X(t) is perfectly unique. 

A good idea to tackle this problem, is to suppose that some of the OD couples do 
not accept any traffic. That is, that they remain equal to zero. The goal is to reduce 
the number of positive OD couples so as to get a system whose routing matrix A(t), is 
square or not too far. 

The first approach is to suppose that some pre-determined OD couples are ex- 
culded. 

To begin, we may partition, a little arbitrarily, the set of the OD couples into two 
parts. In the first one, lie the OD couples which remain always equal to zero. And, in 
the second part, we will suppose that there is some traffic flowing through these OD 
couples. 

We need to generate a sample Z. The components of Z are independent of each 
other and generated from a Bernoulli density whose parameter is pre-determined. Then, 
if Z(i) = 0, the OD couple number i, do not accept any traffic. 



The time required to perform this simulation is of about 1 minute. 

A second point of vue should be to suppose that we know at the beginning that only 
K OD couples, K > n, are positive. So, we need to modify our CE algorithm. We 
need to introduce a matrix, Z(t): 

zP ■■■ z[ N) ] 

7 (1) 7 {N) 
^2 " ' ^2 

(22) 

An ' ' ' A n 

To be more explicit, the i th row of Z is generated from B(pij(t)), j £ {0,1}, Vz 6 
{l,...,n} . In fact, each row is generated independently from a Bernoulli density 
whose parameter is specific, conditional upon the fact that X}j=i -^-j = K, i £ 
{1, N} . K, is a fixed number. It may be as we have already stated, a certain prop- 
portion of OD couples, but it may also take into account some other constraints. 



Z(t) = 
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Figure 6: Introduction of constraints on the OD couples: the zero OD couples are 
pre-determined. Identifiability of X(t) with | of the OD couples used. 

The first idea to deal with such a constraint is to generate a random vector x[ l \ X^ . 
Each component are drawn independently from a Bernoulli density. The sample is ac- 
cepted iff, Y?j=i ■^■j^ — K, i E {1,...,N} . However, when n becomes higher 
than 10, it takes a prohibitive time! In fact, the best solution is to generate independant 
Bernoulli random variables from B{p\(t)), B(p2(t)), respectively, until K unities 
or n — K zeros are generated. Then, the remaining elements are put equal to zero or 
one, respectively. However, the updating formula for the parameters of the Bernoulli 
densities remain exactly of the form: 

Pi,i(t) = — jj x — -,Vi 6 (23) 

In fact, now, X(t) and Z(t) are closely linked. Indeed, if z\ j) = then, x[ j) = 
0, which means than there is no traffic on the OD couple number i for the j th -samp\e. 
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Figure 7: Simulation, the number of zero OD couples is pre-determined, K = |. 



The time required to perform this simulation is of about 2 minutes. 



5 Dynamic estimation 

We have previously determine an estimator of the amounts of traffic flowing through 
the OD couples at the specific instant t. We should ask ourselves whether it is possible 
to determine the trajectories associated to the vector X(t). Particle filtering appear to 
us to be an interesting approach. 
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Figure 8: Identifiability of X(t), K = §. 

5.1 Presentation of Particle filtering 

Particle filtering is a well-known technique based on sequential Monte-Carlo approach. 
It is a technique for implementing a recursive bayesian filter by Monte-Carlo simula- 
tions. The key ideea is to represent the required posterior density function by a set 
of random samples with associated weights and to compute estimates based on these 
samples and weights. 

We will generate a random measure {Cg. fc , Wfc}i=i,...,M that characterises the pos- 
terior pdf p(Co-.k | Xi-k). {CQ. fc ,i = 0, ...,M} is a set of vectors with associ- 
ated weights {w l k ,i — 1, ...,M} (the weights are themselves vectors of weights). 
C*o:fc = {Cj,j ~ 0, is the set of all states up to time t. The weights are nor- 
malised such that, 

M 

= i, Vj e {l •>]■ 

i=l 
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Then, the posterior density at t can be approximated as: 



M 

p(C Q :k | X 1:k ) « J2 W k 6 ( C °* - C 0:k) ■ (24) 

The weights are chosen using the principle of Importance Sampling. Let C % ~ 
q(.), i = l, M be samples generated from a proposal q(.), called Importance sam- 
pling density. By successive approximations, it is shown in [7] that the weights are 
recursively obtained by the following formula: 

wl oc wl i - . (25) 

q(Cl\Cl_ v X k ) 

It can be shown that as M — > oo, the approximation (|24|) approaches the true 
posterior density p(C k \ Xi-.k). 

However, there is a major drawback to use particle filtering techniques. Indeed, a 
common problem is the degeneracy problem. After a few iterations, all but one particle 
will have negligable weight. It has been shown that the variance of the importance 
weights can only increase over time, and thus it is impossible to avoid the degeneracy 
phenomenon. This degeneracy implies that a large computational effort is devoted to 
updating particles whose contribution to the approximation to p(Ck\Xi : k) is almost 
zero. 

Consequently, we have to use resampling mechanisms. The basic idea behind resam- 
pling is to eliminate particles which have small weights and to concentrate on particles 
with large weights. 



5.2 State model and observation equation 

The traffic flow will be modelled as a stochastic hybrid system with discrete states. The 
observation equation is rather simple to get. Indeed, we have: 

X k = Z(C(k)) . (26) 

Where, S is a quite complex function which represent the whole algorithm. 
The difficulty now, is to build a state model. Suppose the flow can be decomposed in 
small entities (for example packets). We note: {Qi> t k\l' S {arcs of the network}, i — ► 
I'}, the number of packets going out of the arc i, during the time interval [ifc,tfc+i[. 
{Qi,k\l G {arcs of the network}, I — > i}, is the number of packets arriving on the arc 
ion [t k ,t k+1 [. 

The conservation of the flow lets us write: 

Yi(k + 1) = Yi(k) + ® l < k - H 0''.* • (27) 

{arcs >z} {arcs V\i—>1'} 
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In fact, 



Qi,k = min(S %tk ; Ri,k+i) ■ (28) 

is called sending function. It expresses how many among the Yi(k) packets in 
the arc i at k are at a distance less than a fixed boundary called f3. Suppose the interac- 
tion between the packets is negligible and their location is uniformly distributed over 
the arc. Si t k is then a random binomial variable with Yi(k) drawings, with probability 
of success arci ^ ngth , or an approximation, since we don't know exactly the length of 
the arc number i. 

The receiving function is defined by: 

R,, k+ i = i Y i max ( k ) + Qhk+i - Yi(k)} ■ (29) 

{arcs l\i — >l} 

The sending function is calculated at first by forward recursion, and we substitute 
Qi.k = Si t k in (|27|) . With this first guess of the amount of traffic in arc i, at time tk+i, 
a first guess of the receiving function can be computed, recursively. Finally, we get: 



Ci(k + 1) = F(Yi(k + l)) = 90ri(k),W(k + l)),Vi e {l,-.,n}. (30) 
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Figure 9: Temporal evolution of the distribution of the traffic on four OD couples. 



6 Conclusion 

We could observe that the performance of the CE method is directly proportionnal to 
the ratio: 

kpt nf nndpsl 

(31) 



| set of OD couples | 

For rather small networks, eg. networks composed of at most 30 nodes, the CE 
method works pretty good and suprisingly fastly. What's more, it is possible to add 
some constraints which could guarantee the identifiability of the vector containing the 
amounts of traffic on the OD couples. At the end of the estimation part, we get estima- 
tors of OD volumes and implicitly, of the routing matrix. In fact, this application is a 
great proof of the simplicity and versatility of the CE method. 

However, some points remain difficult to tackle. For example, when the ratio becomes 
larger than 19, the CE method performs rather poorly. Furthermore, R. Rubinstein 



17 



OD traffic on the network. 1=1 



OD traffic on the network, 1=34 



4- 

3 

2- 
1 



- 







N 

M 
i i 

! \ 



i 

ii 
n 
1 1 

1 1 



1 i;; 

' i.' s 
I I '■■ 



10 15 



20 



1.5 



0.5 



i I 



i I 



1 i \ 
\ i '■ 



n 

■ i 



10 



15 



20 



0.B ■ 



0.6 



0.4 



0.2 ■ 



OD traffic on the network. t=G9 



ii 
ii 
1 1 

I 1 



i ', I 



■■. i 



i i 



10 15 



20 



OD traffic on the network, t=9& 



If <1 i 



10 15 



20 



Figure 10: Spatial evolution of the distribution of the traffic at four different instants. 

recommand that the sample size of the CE algorithm should be of the form: 

N = k n, 5<k< 10. 

Suppose for example, that that we have to deal with a network of 50 nodes. Then, at 
each step of the algorithm we will have to generate a sample of 2450 * 20000 vectors. 
Which is completly impossible due to the limited capacities of our computers. But, it is 
certainly possible to improve the algorithm so as to adapt dynamically the sample size 
to solve this problem. Nevertheless, the question remains open. Fortunately, in every 
network, some specific constraints need to be taken into account. These constraints aim 
at decreasing the number of unknown parameters. The idea to impose that some OD 
couples remain equal to zero is an approach, but there are many others. For example, 
we may want to maximize the global entropy, or other common criteria. 
Particle Filtering is an efficient and subtle technique to dynamically predict the evolu- 
tion of the distribution of the traffic on the OD couples for rather small networks. 
The approaches we use are rather simple to put in application. Nevertheless, they rely 
on many small parameters which are quite difficult to optimize. Furthermore, the size 
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Distribution of the OD traffic. t=1 



Distribution of the OD traffic. t=3S 
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